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ABSTRACT 

Results of a simulation of synchrotron-self Compton (SSC) emission from a 
rotation-powered pulsar are presented. The radiating particles are assumed to 
be both accelerated primary electrons and a spectrum of electron-positron pairs 
produced in cascades near the polar cap. They follow trajectories in a slot gap 
using 3D force-free magnetic field geometry, gaining pitch angles through reso- 
nant cyclotron absorption of radio photons, radiating and scattering synchrotron 
emission at high altitudes out to and beyond the light cylinder. Full angular 
dependence of the synchrotron photon density is simulated in the scattering and 
all processes are treated in the inertial observer frame. Spectra for the Crab and 
Vela pulsars as well as two energetic millisecond pulsars, B1821-24 and B1937+21 
are simulated using this model. The simulation of the Crab pulsar radiation can 
reproduce both the flux level and the shape of the observed optical to hard X-ray 
emission assuming a pair multiplicity of M + = 3 x 10 5 , as well as the very-high- 
energy emission above 50 GeV detected by MAGIC and VERITAS, with both 
the synchrotron and SSC components reflecting the shape of the pair spectrum. 
Simulations of Vela, B1821— 24 and B1937+21, for M + up to 10 5 , do not produce 
pair SSC emission that is detectable by current telescopes, indicating that only 
Crab-like pulsars produce significant SSC components. The pair synchrotron 
emission matches the observed X-ray spectrum of the millisecond pulsars and 
the predicted peak of this emission at 1 - 10 MeV would be detectable with 
planned Compton telescopes. 


1. Introduction 

The recent detection of pulsed emission from the Crab and Vela pulsars at energies 
above 50 GeV by ground-based air-Cherenkov telescopes marks a turning point both in 
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technological achievement of the detectors and in pulsar acceleration and radiation modeling. 
The MAGIC telescope first detected pulsed emission from the Crab above 25 GeV (Alin et 
ah 2008), followed by pulsed detections above 100 GeV by VERITAS (Alin et al. 2011) and 
up to 400 GeV with MAGIC (Aleksic et ah 2012). The averaged and phase- resolved spectra 
above 25 GeV are consistent with a smooth broken power-law extension of the spectrum 
measured by the Fermi Large Area Telescope (LAT) above 100 MeV, thus ruling out an 
exponentially cutoff power-law spectrum predicted by curvature radiation (CR) models for 
GeV emission (e.g. Romani 1996, Harding et ah 2008 [HSDF08]). Very recently, pulsed 
emission was detected from the Vela pulsar above 30 GeV with the HESS (Stegmann et ah 
2014) telescope and above 50 GeV with the Fermi - LAT (Leung et ah 2014). In this case 
it is not clear whether another emission component in addition to CR is required. Most 
recently, VERITAS has reported a non-detection of the Geminga pulsar above 100 GeV 
(Alin et ah 2015). McCann (2015) performed a stacking analysis of 115 Fermi pulsars from 
the Fermi - LAT 2nd Pulsar Catalog (Abdo et ah 2013), both young and MSPs, and found 
no significant emission at more than 7% Crab level above 50 GeV. 

For many years, CR was the favored mechanism for the high energy (> 100 MeV) pulsed 
emission in a number of different models for pulsar emission, and VHE emission similar to 
what has been detected from the Grab was never predicted. In models where acceleration 
and emission occur above the pulsar polar caps (PC models, Daugherty & Harding 1996), 
the predicted spectrum of CR from primary electrons has a very sharp “super-exponential” 
cutoff due to attenuation by the magnetic pair production process that occurs in the very 
strong magnetic fields near the neutron star surface. Measurement of the spectral cutoff of 
the Vela pulsar by Fermi (Abdo et ah 2009) has shown that a super-exponential cutoff can 
be eliminated with high significance, ruling out emission from near the neutron star PCs. 
Models for emission originating in the outer magnetosphere are therefore favored, since many 
of these predict curvature radiation spectra from primary particles with exponential spectral 
cutoffs at energies of a few GeV that agree well with those measured for most pulsars by 
Fermi. Such models include the outer gap (OG), where acceleration and emission occur 
in vacuum gaps that form above the null charge surface (e.g. Romani 1966, Cheng et ah 
2000, Hirotani 2006) and the slot gap (SG), where acceleration and emission occur in a 
narrow region along the last open field lines from the neutron star surface to near the light 
cylinder (Muslimov & Harding 2004 [MH04], HSDF08). More recently, models postulating 
high-energy emission originating outside the light cylinder, near or in the current sheet that 
forms in the spin equatorial plane have been developed. Several of these models assume 
that synchrotron radiation (SR) from particles accelerated by magnetic reconnection in the 
current sheet produces the GeV emission (e.g. Petri 2012, Uzdensky & Spitkovsky 2014), 
while others assume that CR by particles accelerated by induced electric fields in dissipative 
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magnetospheres produces the GeV emission (Kalapotharakos et al. 2014, Brambilla et al. 
2015). 

The Fermi LAT has to date detected more than 160 7-ray pulsars and all of these with 
statistics adequate for spectral fitting of a power law with exponential cutoff 

dN(E)/dE = A 1-7 exp [-(E/E c ) b ], (1) 

show cutoffs, E c , in the range 1-7 GeV (Abdo et al. 2013), assuming b — 1. For a subset of 
pulsars having higher counts, fits to the above expression with b free can be performed and 
a number of these are best fit with b < 1, including the Grab and Vela pulsars, indicating 
a more gradual cutoff than pure exponential. This finding raises the question of whether 
the more gradual high-energy cutoffs are due to a blending of CR spectra from a number of 
different rotation phases (Abdo et al. 2010) and/or locations in the magnetosphere (Leung et 
al 2014), or additional emission components. It has been suggested that the very- high-energy 
(VHE) emission from the Crab is inverse-Compton emission (Lyutikov et al. 2012), since it 
is difficult to produce photons above 100 GeV with CR. Lyutikov (2013) modeled the Crab 
X-ray and 7-ray spectra as cyclotron-self Compton emission from electron-positron pairs in 
the outer magnetosphere, fitting the data to determine the spectrum and multiplicity of 
the pair spectrum as well as the location of emission. An synchrotron-self Compton (SSC) 
model for the Crab VHE emission from the outer gap (OG) was suggested by Aleksic et 
al. (2011) and Du et al. (2012) modeled the observed Crab emission from the annular gap. 
Alternatively, it was proposed that particles accelerated by reconnection in the current sheet 
reach temperatures of 10 GeV, if equipartition is assumed, and their SR can extend to 100 
GeV in the observer frame after Doppler boosting (Uzdensky & Spitkosky 2014, Iwona & 
Petri 2015). 

In this paper, we apply a pair cascade simulation coupled to an outer magnetosphere 
radiation model to predict the broadband X-ray to VHE 7-ray spectra of several classes of 
pulsar. In this model, the pairs are produced in cascades above the PCs initiated by primary 
accelerated electrons and sustained by magnetic pair creation (e.g. Daugherty & Harding 
1982). The pairs lose their momentum perpendicular to the magnetic field to SR near the 
PCs and then stream into the outer magnetosphere, where they absorb radio photons in the 
cyclotron resonance to regain pitch angles and emit SR (HSDF08). The SR provides soft 
photons for inverse-Compton scattering by these same particles to boost their energies to the 
7-ray range. We also include CR radiation of primaries accelerated by electric fields in a SG 
geometry, as well as their SR components due to cyclotron resonant absorption. The radio 
emission is simulated using an empirical cone and core beam model. Since the magnetic held 
structure in the outer magnetosphere becomes important for these models, we use a force-free 
magnetosphere configuration to define the magnetic held, a significant improvement over the 
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retarded vacuum dipole field that we had used previously (HSDF08). Use of the force-free 
magnetic field structure allows us to follow particles and radiation arbitrarily close to and 
even beyond the light cylinder, which should not be a physical limitation of emission models. 
Although force-free magnetospheres assume that the electric field parallel to the magnetic 
field, E || , vanishes, they well approximate the field structure of energetic pulsars where 
departures from force-free conditions exist only in narrow accelerator gaps. We therefore 
use the force-free field structure but assume acceleration of primary particles in a narrow 
SG. Using these simulations, we can address the question of how many pulsars should have 
detectable SSC emission and how high in energy this emission extends. 


2. 3D Magnetosphere Geometry 

Since radiation of relativistic particles is beamed along their direction of momentum, and 
particle trajectories in pulsar magnetospheres follow closely (but not exactly) the magnetic 
field lines, the structure of the magnetic field is a critical component to emission models. A 
static dipole is not a good approximation to a pulsar magnetic field, even near the neutron 
star, since the rotation causes a sweepback of the field lines near and beyond the light 
cylinder. The vacuum retarded dipole solution (Deutsch 1955), adopted in many radiation 
models (e.g. Romani & Yadigaroglu 1995, Cheng et al. 2000, HSDF08), incorporates the 
sweepback due to retardation of the field and produces a distortion and shift of the open 
field volume (including the PC) toward its trailing edge (Dyks & Harding 2004). Numerical 
solutions of force-free magnetospheres (Spitkovsky 2006, Timokhin 2006), where the plasma 
density is assumed to be high enough to screen Ey so that the ideal MHD condition E • B = 0 
is enforced, show a larger degree of magnetic field sweepback, as well as straightening of the 
poloidal field lines due to currents. Recent simulations of dissipative pulsar magnetospheres 
(Kalapotharakos et al. 2012, Li et al. 2012), that drop the ideal MHD condition and allow 
finite Em, show that the structure of the magnetic field stays close to force-free for the high 
plasma conductivity needed to match the y-ray light curves (Kalapotharakos et al. 2014, 
hereafter KHK14) and spectra (Brambilla et al. 2015) of young pulsars. 

We will therefore adopt the force-free magnetosphere structure to model the particle tra- 
jectories and radiation, and all calculations are done in the non-rotating, inertial observer’s 
frame (IOF). The pairs are assumed to be produced at the PCs and to screen the Ey at 
higher altitudes, except for a narrow gap along the last open field lines (the SG). They will 
therefore not experience any acceleration in the calculation of their radiation in the global 
magnetosphere. The primary electrons are assumed to be accelerated only in the SG, with a 
constant Em that we take as a free parameter. Both pairs and primary electrons are injected 
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at the neutron star surface at footprints of open field lines that are determined by “open vol- 
ume coordinates” (Dyks et al. 2004 [DHR04]). These coordinates map the open held region 
bounded by held lines that close within the light cylinder, i?LC = c/O, where hi is the pulsar 
angular rotation frequency. The magnetic held vector is determined at each point by 3D 
interpolation of a table of values read-in from a numerical force-free magnetosphere solution 
(Kalapotharakos et ah 2012). This solution has a resolution of 0.02 i?LC an d extends to a 
minimum radius of 0.2i?LC> which is the surface of a pulsar with period 0.8 ms. To extend 
this solution to the surface of pulsars with longer periods, we join the numerical solution to 
a retarded vacuum dipole solution over the range 0.2 — 0.4i?LC using a ramp function. 

The code hrst computes the PC rim at the neutron star surface by performing 4th order 
Runge-Kutta integrations along the held lines of the numerical solution to determine the open 
held footpoints. This is done by iteration, choosing an initial value of magnetic polar angle 
9-i at a number of different values of magnetic azimuth <fi. If the held line with that footpoint 
does/does not close within Rlc, a smaller /larger value of 9 is chosen until the closure point is 
at i?Lc within a given tolerance. We define “open volume coordinates” (r ovc ,/ ovc ) to identify 
footpoints at the neutron star surface of held lines along which the particles are traced. The 
“radial” coordinate r ovc , equivalent to magnetic polar angle, is equal to 0 at the magnetic 
pole and 1 at the PC rim. Lines of constant r ovc define a set of concentric deformed rings 
that hll the PC surface (see Figure 2 of DHR04). The “azimuthal” coordinate / ovc measures 
the arclcngth along each distorted ring with hxed r ovc . The l ovc increase in the direction of 
magnetic azimuth, <f> pc , in a counterclockwise direction around the PC, starting from l ovc = 0 
at (j) pc = 0, defined to be at the magnetic meridian (i.e. the line between the magnetic and 
spin axes). With this definition, <f> pc = 37t/ 2 is on the leading side of the PC and <p pc = 7t/ 2 
is on the trailing side (Note that <f> pc = 0 + n). 

To model an injection of particles that is uniform over the PC, we use a set of rings 
between r™ r n and with equal spacing d ovc . These footpoints are then spaced uniformly in 
each of the rings with iVj = A az i m 360 equal divisions. Since each ring has the same number of 
footpoints, and the rings have varying circumferences, the contribution from different rings 
must be weighted by hing/Cim, where / ring is the length of a particular ring and / rim is the 
total PC rim length. To determine the trajectories of particles in the IOF, we require the 
total particle velocity to be the sum of a drift component and a component parallel to the 
magnetic held (KHK14), 

( ExB e B\ 
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where (Gruzinov 2008, Li et al. 2012), 


E 2 0 = Bl - B 2 + E 2 

(3) 

Bl = (B 2 - E 2 + y/ (B 2 — E 2 ) 2 + 4(B • E) 2 )/2 

(4) 

B 0 = sign(B • E)JbI. 

(5) 


By requiring that v ~ c and that the motion is outward, we can solve for the scalar quantity 
/ at each point along the trajectory. In the special case of a force-free magnetosphere, v 
stays parallel to B in the corotating frame. In the IOF, the step size along the particle 
trajectory is dd = ds/ f, where ds is the step size along the rotating magnetic held line (see 
Appendix C of Bai & Spitkovsky 2010). In treating the particle trajectories in the IOF with 
force-free held geometry, it is possible to extend the trajectories and radiation outside the 
light cylinder. Although the particles are following the held, they slide forward along the 
held lines as the held lines sweep back, following a nearly radial path in the IOF at r > i?Lc- 


3. Pair Cascade Spectra 

The process of electron-positron pair production in pulsar magnetospheres is thought 
to be critical for generating charges for the magnetosphere and for the pulsar wind nebula, 
as well as plasma necessary for coherent radio emission. The pairs are produced in elec- 
tromagnetic cascades above the PCs (Daugherty & Harding 1982) or in OGs (Cheng et al. 
1986). Here, we assume that the pairs that radiate SR and SSC in the outer magnetosphere 
originate from PC cascades, that are initiated in the strong electric helds near the neutron 
star surface by acceleration primary electrons. We have calculated the spectra of pairs using 
a code that simulates a steady (non-time-dependent) electromagnetic cascade above the pul- 
sar PC (Harding & Muslimov 2011 [HM11]). Primary particles are accelerated by an electric 
held induced by rotation of the magnetic held, derived assuming space-charge-limited how 
(i.e., free emission of particles from the neutron star surface) (Arons & Scharlemann 1979), 
and emit curvature radiation. The highest energy curvature photons are absorbed by mag- 
netic pair attenuation (Erber 1966, Daugherty & Harding 1983), producing a spectrum of 
first-generation electron-positron pairs. The pairs are born in excited Landau states (with 
non-zero pitch angles) and radiate synchrotron photons that spawn further generations of 
pairs. Since the total cascade multiplicity M + (average number of pairs produced by each 
primary particle) depends on the pulsar period P and surface magnetic held strength B s , 
younger pulsars produce high M + cascades but older pulsars, with low magnetic helds and 
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long periods, have much lower M + . This results in a pair death line in P-P which, in the 
case of a dipole field, predicts that the older half of the pulsar population cannot produce 
significant pair multiplicity. 


However, the sweepback of magnetic field lines near the light cylinder (due to real and 
displacement currents), as well as possible field distortions near the neutron star will cause 
the magnetic PCs to be offset from the dipole axis (e.g. Dyks & Harding 2004). HM11 
introduced a perturbed dipole field structure with non-axisymmetry and showed that the 
resulting offset PC enhances pair multiplicity by increasing the accelerating electric field. 
For the pair cascade simulations, we will adopt this magnetic field structure which has two 
configurations for the dipole offset in which the magnetic field is symmetric or asymmetric 
with respect to the dipole axis. Modeling of the thermal X-ray light curves of several MSPs 
requires an offset of the PC hotspot that is symmetric (Bogdanov et al 2013). We adopt the 
symmetric distorted field in simulating the pair spectra of the MSPs in this paper. In this 
case, the magnetic field in spherical polar coordinates (77, 9, 0) is 


„ B s 
B « — 

T) a 


v cos 6 + - 9 (1 + a) sin 9 — (fie sin 9 cos 9 sin (0 — 0 O ) 


( 6 ) 


where r/ = r/R is a dimensionless radial coordinate in units of the stellar radius R , B s is the 
surface field strength at the magnetic pole, a = e cos (0 — 0o) is a parameter defining the 
azimuthal distortion of the polar field lines from a dipole, and 0o is the magnetic azimuthal 
angle that defines the plane of the PC offset. HM11 derive the parallel component of the 
electric field, E», using this field structure and we have used the En of Eqn (11) of HM11 
that corresponds to a symmetric offset. We use this E» to accelerate the electrons from 
different starting points over the PC surface and simulate the pair cascades over the whole 
PCfor a range of P, P (or equivalently, B s ), and offset parameter e. We assume a non-zero 
offset only for the MSPs. Pair spectra for the parameters of the Crab (e = 0) , Vela (e = 0) 
and MSPs PSR B1821-24 (e = 0.6, 0 O = 3vr/2) PSR B1937+21 (e = 0.6, 0 O = 3 tt/ 2) are 
shown in Figure [l] The Crab pair spectrum extends over five decades of pair energy, from 
7 + = 20 to 7+ = 10 6 , while the Vela pair spectrum has a smaller range, cutting off around 
7 + = 4 x 10 5 . The difference is due to the shorter period and larger PC size of the Crab, 
which gives smaller radii of curvature of the field lines, so that CR photons have higher 
energy and gain large angles to the field, and the higher magnetic field strength. The pair 
spectra of the MSPs start at energies about two decade higher, 7+ ~ 2 x 10 3 and extend to 
higher energy, 7+ ~ 10'. This striking difference in pair spectra of young and MSPs comes 
from their large difference (four orders of magnitude) in surface field strengths. As described 
by HM11, the photons must have higher energies to produce pairs in lower magnetic fields, 
which raises the minimum pair energy, but are also able to escape at higher energies, raising 
the high energy cutoff in the pair spectrum. The non-zero PC offset decreases the lower 


cutoff of the pair spectrum relative to the pair spectrum with zero offset by a fraction of a 
decade. 

Simulations of PC cascades including time dependence of the electromagnetic fields 
and plasma production (Timokhin 2010, Timokhin & Arons 2013) find that to satisfy the 
charge and current density of the global magnetosphere models the pair cascades are non- 
steady. The accelerating electric held is time-dependent, in the general case where the 
magnetospheric current density J is not very close to the Goldreich- Julian value, Jqj = Pgjc, 
with cycles of pair creation followed by complete screening of the electric held. The timescale 
of these pair cascade cycles is very short (of order the light crossing time across the gap). 
Such non-steady pair cascades can produce pair multiplicities that are much higher than 
the steady cascades. For young pulsars such as the Crab and Vela, the time-averaged pair 
multiplicities are about several times 10 5 (Timokhin & Harding 2015), compared to ~ 2 x 10 4 
for the Crab and ~ 6 x 10 3 for Vela steady pair cascades. Time-dependent pair cascade results 
are not yet available for MSPs, where 2D simulations are required. Although it is beyond 
the scope of this paper to simulate time-dependent pair spectra (especially for millisecond 
pulsars [MSPs]), we will allow for pair multiplicities higher than that of steady cascades in 
our SSC simulations. 


4. Simulation of Emission 

We model the high-energy radiation over the entire spectrum from optical to VHE 7- 
ray wavelengths. Emission from both primary particles and pairs are simulated as they 
move along their trajectories with step size AE, starting at the neutron star surface to a 
maximum radius r max . Primary electrons are assumed to undergo continuous acceleration 
with a constant electric held, d'f/dE = R acc in a narrow gap along the boundary of the open 
held region, between r ovc = 0.95 — 0.99. The electric held in the high-altitude gap is not the 
same held from HM11 that is used for the pair cascades, since the HM11 held is only valid for 
low altitudes near the PC. The high-altitude gap forms at the outer edge of the PC, where 
the interior E\\ decreases to values too low for pair cascades to develop. Although analytic 
solutions exist for both the low altitude SG (Muslimov & Harding 2003) and its extension to 
high altitudes (MH04), we have chosen not the use the solution of MH04 for several reasons. 
First, the gamma-ray luminosity modeled using the SG E\\ solution of MH04 with widths 
narrow enough to produce the observed light curves falls short of the observed luminosities 
by about a factor of 10 for young pulsars. As noted by Pierbattista et al. (2012), the E\\ 
solution for the original SG and also the OG do not provide enough gamma-ray luminosity 
to match the observed Fermi pulsars. The E n of HM2011 on the offset side of the PC is 



9 


larger and may produce a large enough Eu at high altitude to produce luminosities matching 
the data. However, it is very difficult to derive the Eu extension for the offset PCs and 
this has not yet been done. Also, with the advent of dissipative magnetosphere models 
we will soon be able to have more self-consistent E\\ throughout the magnetosphere, albeit 
numerical. Until these are available, we decided it is better to treat Eu at high altitudes as a 
free parameter. We have therefore simply assumed a constant Eu for the high-altitude gap, 
similar to the behavior of the high-altitude symmetric SG MH04. 


The electron-positron pairs are assumed to experience no acceleration as they flow along 
their trajectories in the screened region inside the gap, between r ovc = 0.91 — 0.95. CR, SR 
and SSC radiation are simulated for all particles in the same way, as described below. The 
primary electrons are injected at the neutron star surface with very low Lorentz factors, 
7 = 2, while the pairs are injected with the pair cascade spectra shown in Figure [lj All 
particles are assumed to initially have momentum only parallel to the magnetic field, since 
the pairs lose all their perpendicular momentum to SR very near the neutron star surface in 
the pair cascade. Single primary electrons are injected at the start of each primary trajectory 
at the neutron star surface and their radiation is normalized to the primary flux 


n p = n GJ ckR 2 9 2 pc [(C ») 2 - (C“ ) 2 ] (7) 

where uqj = B 0 Q/2Tcec is the Goldreich- Julian density and 9 P c = (VtR/c) 1 ^ 2 is the PC 
half-angle. The primary particles from the interior of the PC are on field lines that do not 
thread the high- altitude accelerator gap, but are on field lines where the accelerating field is 
screened above a pair formation front by pair cascades at low altitude. Since these primaries 
are not accelerated above the pair front, they lose most of their energy quickly and do not 
contribute much to emission at high altitudes. The emission from the PC pair cascades is 
all emitted at low altitudes, below 6-7 stellar radii, and will only be visible when an observer 
viewing angle is very close to the magnetic axis. Since most observer angles will not see the 
PC emission, we therefore can neglect their contribution in the present study. 


To simulate radiation from a spectrum of pair energies, shown in Figure [lj the pair 
spectrum is divided into logarithmic energy intervals, A 7 + , between 7 ™ m and and a 
single pair member (electron or positron) at each energy, 7 +, is injected at the base of each 
pair trajectory. Their radiation is normalized to the flux of pairs in each energy interval, 


w+(7+) 


2M + n + (7 + )A7+hp 


~.max 

f^in n+(7+)d7+ 


( 8 ) 


and M + is the pair multiplicity. We have kept M + as an independent quantity rather than 
derive it from the integral of the pair spectrum so that it can be treated as a free parameter 
in the simulation. 
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The direction of photon emission, rj em , is along (3 — v/c, the direction of particle motion 
in the IOF, with photon emission angles 


Hem 


Pz, 


(pern — cttcLR 



( 9 ) 


As the particle moves a distance d£, the change in phase is A0 rot = Qd£/c = Qdt. The 
emitted photons are accumulated in sky maps of observer angle £ 0 b s = acos (/x em ) versus 
phase 4> 0 hsi both with respect to the pulsar spin axis (the z axis on the IOF Cartesian grid). 
To determine the observed phase, we apply the rotation and add the time delay correction, 

</>obs = 0em - A0 rot + rem 776771 (10) 

me 

where r em is the radius of photon emission. The quantities in the sky maps are the emitted 
photon fluxes, normalized using Eqn (JtJ) for primaries and Eqn (J8]) for pairs, divided by the 
solid angle of that sky map bin, = sin£ o b s A£ o b s A0 o b s . The phase-averaged observed flux 
at a viewing angle £ 0 b s is then obtained by summing the fluxes in </> 0 b s and dividing by 2tc d 2 , 
where d is the distance to the source. 


4.1. Curvature Radiation 


The energy spectrum of curvature radiation from a single electron with Lorentz factor 


7 is 


A [ cr(£) = ^3 — 7 k ( — \ 

C \£cr / 


where e is the emitted photon energy in units of me 2 and 

3 c o 
= X— 7 , 

2 pc 

and the function n(x) is defined as 

poo 

k ( x ) = X / K, 5/3 (x') dx'. 


( 11 ) 


( 12 ) 


(13) 


The radius of curvature p c in the force-free magnetosphere is not that for a pure dipole field, 
but is the radius of curvature of the particle trajectory determined in the IOF by computing 
the inverse of the trajectory curvature using the second derivative at the particle position: 


C 


d 2 v 


d £ 2 


P> 


-l 


(14) 
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The approximate form of the photon spectrum is a power law with an exponential cutoff 


at e 


CT 5 


ds 


a 

(A me) 1 / 3 



e 2/3 exp(-e/£ or ). 


(15) 


4.2. Synchrotron Radiation 


Since all particles start their trajectories with momentum parallel to the magnetic field, 
they will not radiate any SR until they acquire finite pitch angles. We assume that particles 
gain pitch angles by resonant absorption of radio photons (Shklovsky 1970, Lyubarski & 
Petrova 1998, LP98). In this process, relativistic particles absorb photons that have energies 
at the cyclotron resonant frequency in their rest frame, which results in an excitation of the 
particle to a higher Landau state and an increase in pitch angle. The particle will then emit 
spontaneous cyclotron radiation if its momentum perpendicular to the magnetic field is non- 
relativistic in the frame where the parallel momentum vanishes, or synchrotron radiation if 
the perpendicular momentum is relativistic. The resonant absorption condition is 


B' = 7£ 0 (1 - Pfio) 


(16) 


where 7 is the particle Lorentz factor, £0 is the lab frame energy (in units of me 2 ) of the radio 
photon, /3 = (1 — I/7 2 ) 1 / 2 , B' = B/B c r is the local magnetic field in units of the critical field 
strength B cv = 4.4 x 10 13 G, /i 0 = cos 6 * 0 , and 9 0 is the angle between the photon direction 
and the particle momentum in the lab frame . In pulsar magnetospheres, particles having 
large Lorentz factors will see radio photons at the cyclotron resonance at high altitude above 
the neutron star surface, when the local magnetic field has dropped low enough to satisfy 
the resonant condition, (from equation |T6] ) , 


7 R 


= 2.8 x 10 5 


B 8 

£ o , ghz (1 — Pno) 


( 17 ) 


where /x 0 is the incident absorption angle. 


The particle initially undergoes absorption in low Landau states, where the cyclotron 
emission rate is well below the absorption rate. Therefore, the particle Landau state (and 
pitch angle) will increase stochastically but continuously until the increase in pitch angle 
through resonant absorption and the decrease in pitch angle by synchrotron emission reaches 
an equilibrium. This balance is not reached until the particle occupies a high Landau state, 
where it will radiate synchrotron emission. 


We use the method of HSDF08, based on the work of LP98, to simulate the resonant 
cyclotron absorption and subsequent synchrotron emission from both primary electrons and 
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pairs. LP98 identified two regimes of particle pitch angle increase in the resonant absorption. 
In the first regime, when 0 <C 9 0 (the particle pitch angle, -0, is less than the radio photon 
incident angle of the, do), the pitch angle is increasing while the momentum is nearly constant. 
In the second regime, when (9q — 0) -C do, the pitch angle is constant as the total momentum 
increases. In the regime where 0 do (equation (2.17) of Petrova (2002)), the mean square 
of the pitch angle increases as, 


(0 2 )=4 / ao{rf)drf, 
dr/R 


(18) 


where r) — r/R, r) R is the radio emission altitude and 

27r 2 e 2 (l - f3no)I 0 {£0^(1 - f3n 0 )y 
ttoM = iW ( B ) ' 


(19) 


Here Jo is the observed intensity of the radio emission, in erg ■ cm~ 2 ■ s _1 ■ Hz 1 and v is 
its spectral index. Thus, the change in perpendicular momentum due to cyclotron resonant 
absorption is 


dp± 

dt 


abs 


o / N 7 ,P± ( dp 

= 2 a 0 {ri) c h — — 

P± P \dt 


abs 


( 20 ) 


where we used the assumption that p± is proportional to the root mean-square of the pitch 
angle, pj_ = p(0 2 ) 1//2 . We also compute the evolution of p(0 2 ) 1//2 rather than computing the 
evolution of the full particle distribution function. 


Combining Eqn (19) and (20), the resonant absorption rate is 


dp± 

dt 


abs 


7 

= £— + 


P± 7 


P± 7 


dt 


abs 


7 < 7 r 


where 


D = 5.7 x 10 9 s" 1 7 ^ 


-^kpc 


^o[mJy] (1 - /3/x 0 ), 


( 21 ) 


( 22 ) 


and we have neglected the (d^/dt) abs term in Eqn (21) since the (d'y/dt) from acceleration 
and from curvature and synchrotron losses are much larger. In Eqn ( [22]) , we take Jo = 
f&o ^radd 2 /A, where Vt ra d ~ A/r 2 is the solid angle of radio emission, with A the cross- 
sectional area and r the radius at absorption, $0 is the measured radio flux (in mJy), and 
d is the source distance (in kpc). If and when the particles satisfy the resonance condition, 


7 < 7 r (where 7 r was defined in Eqn (17)), as they advance along their trajectories, the 


resonant absorption term Eqn (21) will turn on 
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In the regime where 9 0 — i/j <C 9 0 , we have assumed that the pitch angle maintains a 
constant value of ip = 9 0 /2 and that the mean of the particle total momentum is 


P = 



(3 — v)a\ 
M o 


1/(3-!,) 


(23) 


(Petrova 2003), where 


ai {v) 


47r 2 e 2 J , i(l)/ 0 


£ o@o \ 
B'mc ) ’ 


V > V R , 


(24) 


bj = 2e 2 B’ 2 /3h 2 e and J[ is the derivative of the Bessel function. 

The radio emission is modeled as an empirical cone and core beam centered on the 
magnetic pole to determine the flux distribution of radio photons S(9,£r ), given by Eqn (9) 
of HSDF08, where 9 is the magnetic polar angle and £r is the radio photon energy in units 
of me 2 . We assume this core component of the flux is emitted at a radius of 1.8 R and the 
conal component at radius (in units of the stellar radius) 


tk g ~ 40 


P 


10 


-^ss- 1 


0.07 

p0.3--0.26 
r b GHz 


(25) 


(Kijak & Gil 2003), where Eqhz is the radio frequency in GHz. In order to evaluate the radio 
photon intensity <h 0 and the incident absorption angle /i$ at the position of the particle, 
needed in Eqns (22) and (24), we use this radio core/cone beam model. To evaluate <3>o, we 
divide the radio emission beam in ovc coordinates, defined in (|2j into “beamlets” centered on 
tangents to the field lines and having opening angle One set of beamlets is modulated 
by the cone beam at radius tkg (see Eqn (25)) and another set of beamlets is modulated 
by the core beam at radius 1.8 R. The absorption angle /io between each beamlet and a 
particle with which it will interact is determined by computing the travel time of the radio 
photons from their emission points (ay )m ,ybm,Zbm) to the particle position ( x p ,y p ,z p ), taking 
into account the rotation of the field line during the photon transit time. A particle at a 
given location in the outer magnetosphere can absorb only a fraction of beamlet photons. 
Contributions to the cyclotron absorption from all the separate beamlets is summed at each 
particle position. 


We have used this empirical model of radio core and cone emission to model the cyclotron 
resonant absorption, although some of the pulsars we study show phase alignment of their 
main radio and high-energy pulses, indicating that the radio emission comes from high 
altitude near the y-ray emission. However, the radio cone emission altitude from the Kijak 
& Gil (2003) model for the Crab at 400 MHz is about 0.27 Alc and for B1937+21 is about 
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0.55 i?LC; which is at fairly high altitude and would likely form a caustic emission pattern to 
give the phase alignment with the gamma-ray pnlses. Harding et al. (2008) modeled both a 
cone beam at this single altitude and a cone beam with some extension along the field lines 
to higher altitude for the Crab in and the resulting SR was similar for both cases. In the case 
of the Crab and some MSPs like B1821+24 and B1957+20 some radio peaks are in phase 
with the gamma-ray peaks and some are near the phase of one of the magnetic poles (as 
for the Crab precursor), indicating that there is both low altitude cone or core emission and 
high altitude radio emission simultaneously. The cones may also be elliptical, not circular, 
from studies of precessing pulsars (Weisberg & Taylor 2002). We believe our radio model 
is sufficient for the present calculations since it does place the cone beam emission at high 
altitude for the Crab and the MSPs, in the same region with the gamma-ray emission. 

Although the electron-positron pairs from the PC cascades are thought to produce the 
observed radio emission, we have neglected any energy loss of the pairs due to emission of 
the radio beam. This is reasonable since the radio emission power is a small fraction of 
the spin-down power (about < 10 -6 ) compared to the fraction of spin-down power used to 
produce the pairs (about ICC 3 , HM11). 

At each step along its trajectory, the particle radiates an instantaneous synchrotron 
spectrum (Tademaru 1973), 

2 2 / 3 

N S r(£) = ^ I yaR / sin^£" 2/3 £^ /3 exp(-£/£ SK ), (26) 

where sin^ = p±/p, p 2 = y 2 — 1 and e SR = (3/2)y 2 B' sin-0 is the synchrotron critical 
frequency. 


4.3. Synchrotron Self-Compton Radiation 


Calculation of the SSC emission, the inverse- Compton radiation of both pairs and pri- 
mary particles scattering the SR from both pairs and primaries, is carried out in two stages. 
Since the SR photon density is needed to compute the SSC emission, the first stage cal- 
culates the SR from both pairs and primary particles at each step along all trajectories in 


the open field volume according to 14.2 and also accumulates the SR emissivity at each 
location, r em = (x em ,y em , z em ), in Cartesian coordinates, and photon emission direction, 
ijem = {Vem,x, Vem, y , Vem,z) in the IOF. At each position along a particle trajectory, depend- 
ing on whether the particle is a primary electron ( p ) or a pair (+), the SR emissivity is 
incremented by 

N S r{s ) fi p j+ 


Acgj^iS, r em , T ] ern ) 


c Ax em A^/ em A^ 

em 


(27) 
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where Ad is the spatial step along the particle trajectory in the IOF, h p is the flux of primary 
electrons, given by Eqn (JT|) , h + is the flux of pairs from Eqn (J8| summed over pair energies, 
and Ax em , A y errL , A z em are the spatial divisions in the r em array. At the end of this first 
stage, the SR emission and the total SR emissivity, e RR (£,r em ,ri ern ), throughout the open 
magnetosphere of one hemisphere (northern) has been computed. We need the SR emissivity 
from both hemispheres, since particles along trajectories in one hemisphere can scatter SR 
photons produced by particles on trajectories in the other hemisphere, the emissivity from the 
southern hemisphere is computed from the one in the northern hemisphere using reflection 
symmetry: 

^SrYi ^ e mi ^7em) ^SR Y > hnn ^7 'em)- (28) 

The total SR emissivity is then, e SR (£, r em , r] em ) = eg R (e, r em , r/ em ) + ef R (e, r em , r/ em ). 


In the second stage of the simulation, the particles follow their trajectories for a second 
time, radiating SSC using the SR emissivity computed in the first stage. The scattered 
photon distribution from a single particle is, 


dNsscjts) 

de s dt 


J d£l J den ph {e,Si) (1 - d/i ) ' 


(29) 


where 7i p h(e,fl) is the SR photon density, cr(£,h2) is the scattering cross section, and £ and 
e s are the incident and scattered photon energies. At each step along the particle trajectory, 
the SR photon density in all directions is computed at the current position, r IO F, of the 
particle, 


Tlph AlOF ; ^7em) ^ I dx e 


rymax 
' Umin 


dy e 


dz P 


£ sr ( z , ^*eni9 ^7e 


where 

r s = ( x IOF — x emY + (l/lOF ~ Hem Y + (ZlOF ~ Z em ) 2 . 


(30) 


(31) 


In principle, the Klein-Nishina (KN) scattering cross section should be used in Eqn 
(29). In practice, a full integration over all angles using the full KN cross section at each 


particle position is not computationally feasible. We have therefore adopted an approximate 
approach, using the formula of Jones (1968) for the KN photon production rate for a sin- 
gle ultra-relativistic particle with Lorentz factor 7 scattering an isotropic distribution, and 
modulated by the anisotropic photon density and the relative velocity factor: 


Af sAA ^ J. f dn 

a£edt 47 t / 


de n ph (£,rioF,r] em ) (! - Pt*) 


dn K jy(g,e s ) 
(feoCfe 


(32) 


where 
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dn KN (£,e s ) 2 tt r%c 1 ( Yq ) 2 

= ~ 2~ [2glng+ (l + 2g)(l - g) + - - , - _ (1 - g)], 


de.de 




21 + Tg 


(33) 


r = 4£7, 9 = 


E\ e s 

E-t = — 

r(i 7’ 


(34) 


and r 0 is the classical electron radius, // = cos# is the incident photon angle and (3 = 


(1 — I/ 7 2 ) 1 / 2 . In the limit T « 1 and Ei « 1, Eqn (J32J) reduces to the Thompson limit 
result (Blumenthal & Gould 1970). 

To compute the SSC spectral contribution from a particle at position tiof with velocity 


v in the IOF, the integration in Eqn (32) is performed by looping over the whole range 


of incident photon energy £ and directions /1 and 0, relative to the particle direction. We 
assume that the scattered photon direction, fi s and </> s , is the same as the particle direction, 
rje, SO 


, 1 de,y 

Vs = Ve,z, tan (f> s = 

T/e,z 


(35) 


To evaluate the photon density for each incident photon direction we need to know the 
incident photon direction, /q and 0j, in the IOF, 


Vi = VsV + s i n @s sin 6 cos 0 


and 


with 


cos A 0j = 


04 0s T A (pi 

cos 6 — cos 6 b cos 6j 


(36) 

(37) 

(38) 


sin 6 S sin 6i 

The quantities Vi an d 0i define the direction rji in the IOF that is then used to find the SR 
photon density in that direction. 


4.4. Particle Dynamics 


The equations of motion for the Lorentz factor, 7 , and perpendicular momentum, p± 
(in units of me) of a particle as it moves along a field line can be written (Harding, Usov & 
Muslimov 2005), f 

dr [= eE 1 _ Je* 2 _ 2eV ( drA abs _ ( dj_\ SSC 

dt me 3 m 3 c 5 3 p 2 \dt ) \dt J 
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dp± _ 3 c 2e 4 / dp ± (y) \ 

dt 2 r^ >± 3 m 3 c 5 7 \ dt J 


(40) 


In equation (39), the different terms of the right hand side are acceleration, synchrotron 
losses, curvature radiation losses, cyclotron/synchrotron absorption and inverse Compton 
losses. I11 equation (40), the various terms of the right hand side are adiabatic momentum 
change along the field line, synchrotron losses and cyclotron/synchrotron resonant absorp- 
tion. The SSC losses are negligible for p±. 


In the first stage, integrating the particles trajectories to compute the SR and the SR 
emissivity, only SR losses are included since the SSC losses are computed in the second 
stage. However, we find that the SSC losses are much smaller than the SR losses in all cases 
except for the very high altitude part of the pair trajectories for the Vela pulsar (see Figure 
[2] discussed below). For this case, discussed in j|5j the SSC losses are negligible, and the SR 
losses do not change the particle energy significantly since the cyclotron absorption balances 
the SR losses. The spatial step size for the particle trajectories is dynamically adjusted at 
each step to not exceed limits set by radius of curvature of the trajectory (to limit fractional 
change in emission direction to 10%), acceleration rate (to limit fractional gain in primary 
Lorentz factor to 10%), CR and SR loss rates (to limit fractional energy loss to 10%), and 
absorption rate (to limit fractional increase in p± to 50%). 


5. Results 

The radiated spectra of primary particles and pairs have been simulated for the Crab 
pulsar, the Vela pulsar, and two of the most energetic MSPs in the Northern and Southern 
hemispheres, PSR B1937+21 and B1821-24. The pair spectra computed for each of the 
pulsars using the cascade simulation described in j|3] are used to compute the spectra of SR, 
CR and SSC radiation. However, we tried several different pair multiplicities as well as the 
multiplicity from the original cascade simulation. 

The particle trajectories start at the neutron star surface and extend to r max = 2.0Rlc 
in radius or = 1.9Rlc in cylindrical radius, whichever is reached first. The rj^ ax limit 
is set so that the calculation stays inside the grid of the numerical field lines. The region of 
footpoints of injected particles on the PCs is divided into 5 rings between r™™ and r™^ x , each 
of which is divided into 15 azimuthal divisions (A azim = 0.04 divisions per degree), totaling 
75 each of injected primary particles and pairs at each of 34 energies. The code was run on 
75 parallel processors on the Discover cluster at Goddard. For the calculation of the SSC 
spectrum, we used 20 divisions in cos 6 and 12 divisions in </> in the integration over incident 
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photon directions (Eqn [32]). The SR photon emissivity array has 50 logrithmic energy 
bins, from 2 x 10 -6 — 10' MeV, 18 angular divisions in each of three independent Cartesian 
direction cosines, rj ern = (rj em ,x,'nem,y,'nem,z) from —1 to 1, and 9 divisions in each of three 
spatial coordinates r em = (x em , y em , z em ) from — 2.0-Rlc to 2.0Rlc- To assure that there 
was enough angular accuracy in the SR photon emissivity array of the SSC computation, we 
tested the computed SSC flux level with increasing number of angular divisions. It was found 
that the SSC flux decreases with increasing number of divisions but stabilizes at 18 divisions 
in each component of r/ em , for a combination of 18 3 = 5832 elements over 47 T steradians, so 
this number was adopted. 


Figure [2] shows examples of the evolution of the pair perpendicular momentum, p±, 
cyclotron absorption rate, (dp±/ d£) abs , synchrotron loss rate, (d , y/d£) SR , SSC loss rate, 
(dy /d£) ssc , and synchrotron critical energy, £$r , as a function of radial distance along a 
pair trajectory for the Crab, Vela and B1937+21, from integration of Eqns (39) and (40), all 
for pair energies around 7 + ~ 10 5 . The cyclotron absorption starts at the radius of the radio 
cone emission, which is about 0.2i? L c f° r the Crab, about O.IRlc for Vela, and 0.56Rlc f° r 
B1937+21, according to Eqn (25). Before the radio emission altitude is reached, there is 
some transient SR from having to set a small non-zero initial p± = 10 -5 , to avoid overflow 
in evaluation of (dp±/d£) abs in Eqn (21). The fluctuations in the absorption rate as a func- 
tion of radius are due to emission from different beamlets the electron encounters (see 14.2), 
reflecting the numerical resolution of these sub-elements of the radio cone beam. For this 
pair energy, the p_±_ becomes relativistic. In the case of the Crab and B1937+21, the radio 
photons are emitted at high enough altitude and the magnetic field strength near the light 
cylinder, Rlc, is high enough that the radio photons stay in resonance in the particle rest 
frame through the entire particle trajectory. Their Lorentz factors (not plotted) are nearly 
constant since the cyclotron absorption rate comes into equilibrium with the SR loss rate. 
In the case of Vela, the absorption occurs over only a small part of the trajectory above 
the radio emission altitude, above which the particle is 110 longer “in sight” of the relatively 
more narrow radio beam. The p±, (d'y/d£) SR and Esr thus drop with increasing radius. The 
SSC loss rate is always less than the SR loss rate, except in the case of Vela where the SR 
loss rate drops below the SSC loss rate at higher altitude. The SSC loss rate for Vela does 
not decrease much at high altitude. Because the T>lc for Vela is much lower compared to 
the other pulsars, the SR critical frequency is about four order of magnitude lower, so that 
the radiated photon number density increases with altitude. Although the SSC loss rates 
are not included in the trajectory integration, they are never large enough to change the 
particle energy significantly. 


Figure [3] shows examples of the evolution of the primary electron Lorentz factor 7 , 
curvature radiation loss rate, (dy / d£) CR , synchrotron loss rate, (dy /dt) SR and SSC loss rate, 
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(d'y/di ?) ssc , as a function of radial distance along a pair trajectory for the Crab, Vela and 
B1937+21. For the primaries, CR losses are always much larger than SR and SSC losses and 
the electrons reach CR reaction limit in all cases within a few tenths of Rlc- The maximum 
Lorentz factors for the chosen values of acceleration rate are slightly over 10 7 , giving critical 
CR energies of a few GeV. The SR loss rate peaks soon after the start of the cyclotron 
resonant absorption, at which point the particle Lorentz factor drops as the SR loss and 
acceleration gain rates are briefly in equilibrium (Harding et al. 2005). As the SR loss rate 
drops, the particle then returns to equilibrium between acceleration gain and CR loss rates. 


5.1. Crab Pulsar 

For the Crab pulsar, we use the following measured parameters: P = 0.033 s, P = 
4.22 x 10 _13 ss _1 , d = 2 kpc and SHOO = 700 mJy, where d is the distance and SHOO is the 
radio flux at 400 MHz. In addition, we assumed a magnetic inclination angle of a = 45°, a 
constant acceleration rate of R acc = eE\\ /me 2 = 0.2 cm -1 = 6 x 10~ 4 Blc, where Rlc = 6 x 10 5 
G, for the primary particles, and two pair multiplicities, M + = 2 x 10 4 , from a steady-state 
pair cascade, and M + = 3 x 10 5 , from a time-dependent pair cascade. 

Figure [4] shows the observed and modeled spectral energy distribution (SED) of the 
phase-averaged flux of the Crab pulsar from optical to VHE y-ray energies at viewing angle 
( = 60° and for M + = 2 x 10 4 . The different model radiation components, primary CR, 
SR, SSC and pair SR and SSC, are plotted separately. We multiplied the phase-averaged 
pair SR flux by a factor of 11 to match the observed spectral points and multiplied all the 
other components by the same factor. The shape of the pair SR component, dependent 
on the shape of the pair spectrum, and its peak energy, dependent on the magnetic field 
strength near the light cylinder, nicely matches the observed spectrum from optical through 
the BeppoSax hard X-ray measurements. However, it falls short of the COMPTEL low- 
energy y-ray points. The pair SSC component peaks at several GeV but falls well below the 
observed Fermi and VHE spectral points. The primary CR component peaks near a GeV 
but way above the observed spectrum, while the lower primary SR component peaks around 
1 MeV at the level of the observed spectrum. The primary SSC component, peaking sharply 
at a few times 10 6 MeV, has a flux below the scale of the plot. 

Figure [5] shows the observed and modeled SED for the same parameters as the model 
in Figure [4] but for a higher pair multiplicity of M + = 3 x 10 5 . For this case, the pair SR 
spectrum matches the observed spectrum without any renormalization factor. The associated 
pair SSC component now roughly matches the level of the higher-energy Fermi points and 
the MAGIC and VERITAS measurements. For the adopted value of the acceleration rate 
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R acc , the primary CR component fills in at a few hundred MeV. However, the primary SR 
component is too low to account for the observed soft 7-ray emission. The tip of the steeply 
rising, nearly mono-energetic primary SSC spectrum is just visible at a few TeV at the 
bottom right edge of the plot, and the pair CR lies several orders of magnitude below the 
plot lower boundary. 

We tested the addition of a high-energy power law extension to the pair spectrum (shown 
in Figure [5| whose SR spectrum would account for the observed emission in the 1 - 100 MeV 
range. The resulting pair SR and SSC opponents are shown in Figure [5] as dashed lines. 
The SSC spectrum now extends to higher energy with a harder spectrum that exceeds the 
observed MAGIC and VERITAS points. This extension of the pair spectrum, which has no 
physical basis, would thus produce an SSC spectrum that seems to be ruled out by the data. 
This implies that the observed 1-10 MeV emission is not produced by the same particles 
that produce the SSC emission. 


5.2. Vela Pulsar 

For the Vela pulsar, we use the following measured parameters: P = 0.089 s, P = 
1.25 x 10 -13 ss -1 , d = 0.29 kpc and 5400 = 5000 mJy. We assumed a magnetic inclination 
angle of a = 75°, a constant acceleration rate of R acc = eE\\/mc 2 = 2.0 = 0.08 T>lc for 
the primary particles, and two pair multiplicities, M + = 6 x 10 3 , from the steady-state pair 
cascade, and M + = 1 x 10 5 , for a time-dependent pair cascade. 

Figure [6] shows the observed and models SED of the phase-averaged flux of Vela at 
viewing angle Q = 60° for the two different values of M + . We multiplied the phase-averaged 
primary CR flux by a factor of 0.14 to match the level of the observed Fermi spectral points 
and multiplied all the other components by the same factor. The value of R acc had been 
chosen so that the SED peak of the primary CR spectrum matches the peak of the Fermi 
spectrum, but the flux level being too high implies that either the flux of primary particles 
or the distance over which they are radiating is smaller than what we assume. The pair SR 
and SSC components peak at energies of around 1 keV and 1 GeV respectively, somewhat 
lower than the peaks of those components for the Crab. For both values of M + , the SSC 
flux is orders of magnitude lower than the primary CR flux. The primary SSC component 
is well below observable levels. The striking differences between the Vela and Crab SEDs 
come from a number of differences in their intrinsic properties. First, although they have 
similar surface magnetic field strengths, the larger magnetosphere of Vela gives it a much 
lower 7 ?lc = 4 x 10 4 G compared to that of the Crab. This produces the lower peak energies 
of the SR and SSC components and also the lower flux levels because the magnetic field 


21 


drops below the value where the radio frequencies are in the cyclotron resonance in the 
pairs’ rest frames well before the light cylinder, whereas for the Crab, the radio photons are 
in resonance out to and beyond the LC. Second, the radio luminosity of Vela is lower by 
about a factor of 5 compared to the Crab, and is produced at lower altitude relative to the 
LC. The particles therefore see a lower density of radio photons over their trajectories and 
undergo less absorption. Third, the pair spectrum of Vela has a high-energy turnover at a 
lower energy, producing the lower peak energies of SR and SSC spectra. 


5.3. Millisecond Pulsars 

MSPs are promising sources of high non-thermal SR and SSC emission, since they have 
small magnetospheres and consequently a number have large Rlc- In fact, a larger fraction 
of MSPs have observed hard non-thermal X-ray spectra that do non-recycled pulsars. We 
model two of the most energetic MSPs, B1821-24 and B1937+21, that have the highest R L c 
and high levels of non-thermal X-ray emission. B1821-24 (J1824-2452), with P = 3.05 ms, 
P = 1.6 x 10 _18 ss^ 1 and 5400 = 40 mJy, lies in the globular cluster M28, at a distance of 
d = 5.5 kpc. B1937+21 (J1939+2134), with P = 1.6 ms, P = 1.0 x lO^ss" 1 , 5400 = 240 
mJy and d = 3.6 kpc, was the first discovered MSP (Backer et al. 1982). Both are isolated 
MSPs, have B L q = 7.3 x 10 5 G and 10 6 G respectively, exhibit giant radio pulses (Bilous et 
al. 2015, Popov & Stappers 2003) and have aligned high-energy and radio profile components 
similar to the Crab. 

Figure [7] shows the observed and modeled spectra of B1821-24 from soft X-rays through 
VHE y-ray energies for a = 45° and viewing angle ( = 80°. Although there are no strong 
constraints on either a or (, modeling the q-ray light curve only favors a ~ 40° and £ = 85° 
(Johnson et al. 2013). We compute SR and SSC spectra for M + = 10 3 , a value for a pair 
cascade with offset PC of e = 0.6, and also an extreme value of M + = 10 5 , both assuming 
-Race = eE\\/mc 2 = 2.0 = 5 x 10~ 3 Rlc f° r the primary particles. The pair SR component 
peaks between 1 and 10 MeV and the pair SSC component peaks around 50 GeV, almost 
two decades higher than the SR and SSC components in the Crab SED. This is expected, 
given than the pair spectrum of B 182 1-24 is higher in energy by a similar factor relative 
to the Crab pair spectrum, and the B^c values are similar. The pair SR spectrum for 
the M + = 10 3 case matches the slope and level of the X-ray data with a minor shift by a 
factor of 0.8, which has been applied to the other components for this model. The primary 
CR component roughly coincides with the Fermi points, with the primary SR being much 
lower. The pair SSC spectrum lies about four decades below the pair SR and primary CR 
components and is not detectable. The pair SR flux for M + = 10 5 is almost two decades 
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above the observed X-ray flux, although the pair SSC flux might be detectable with Fermi 
and even HESS. Since time-dependent cascade simulations for MSPs are not available at 
present, we cannot say whether such high multiplies are in fact achievable, although a pair 
multiplicity much higher than M + = 10 3 , which already accounts for the observed X-ray 
spectrum, seems to be ruled out. Therefore, for our chosen values of inclination and viewing 
angle, we do not predict a detectable VHE component from B1821-24. 

Figure [8] shows the observed and modeled spectra of B 1937+21 from soft X-rays through 
VHE y-ray energies for a = 75° and viewing angle ( = 70°. Johnson et al. (2014) fit the 
combined y-ray and radio light curves of this MSP using several geometrical emission models. 
For the two-pole caustic geometry, they find a = 88°(+2 — 1)°,C = 88°(+2 — 1)°, and for 
an OG geometry, a = 72° ± 1°,£ = 85° ± 1°. We show the spectra for pair multiplicities 
of M + = 10 3 of a pair cascade with e = 0.6, and also for M + = 10 5 , both assuming 
Race = 2.0 = 3 x ICE 3 Hlc for the primary particles. Applying only the l/(27rd 2 ) factor to 
obtain the phase-averaged flux from the sky map, the pair SR flux for M + = 10 5 is consistent 
with the Chandra X-ray flux. The peak of the pair SR SED is around 5 MeV, very similar 
to the SR SED of B1821-24. The primary CR component is about a factor of 3 below the 
Fermi flux points, and the pair SSC component, peaking at around 100 GeV, is about a 
factor of 2 below CTA sensitivity. The pair SR component for the case of M + = 10 3 is about 
two decades below the observed X-ray flux and the pair SSC emission lies well below the 
plot boundary. The SR fluxes of B1821-24 and B1937+21 for the same pair multiplicity and 
similar are very different since the model inclination angles are quite different. The value of 
a = 75° was chosen to match the best model fits of the y-ray light curve, but a simulation 
for B1937+21 with a = 45° would match the observed X-ray spectrum with a lower pair 
multiplicity. In this case, the SSC component would be much lower. 


6. Discussion 

We have simulated the synchrotron and synchrotron-self Compton emission from a broad 
spectrum of electron-positron pairs in a global force-free magnetosphere. The scattering takes 
place in the extreme Klein-Nishina limit and although the angular dependence of the cross 
section is neglected, the full 3D angular dependence of the SR photon density is treated by 
storing the SR emissivity over the entire magnetosphere from both rotational hemispheres. 

Our simulation of the Crab pulsar radiation reproduces both the flux level and the 
shape of the observed optical to hard X-ray emission assuming a pair multiplicity of M + = 
3 x 10 5 . Such a high multiplicity is about an order of magnitude larger than that produced in 
steady-state PC cascades, but is achievable with the more realistic time-dependent cascades 


(Timokhin & Harding 2015). The predicted SSC flux in this case roughly matches the tail 
of the observed Fermi spectrum as well as the MAGIC and VERITAS detected points. A 
CR component from primary particles is necessary to explain the Fermi spectrum below 
a few GeV. The model SSC spectrum is not a power law, but reflects the shape of the 
pair spectrum, which has a smooth and continuous curvature extending to ~ 1 TeV. An 
SSC component from primary particles appears above a few TeV and may ultimately be 
detectable. 

The Vela pulsar does not produce detectable SSC emission in our simulations, even for 
a pair multiplicity of 10 5 . The predicted SSC emission has both a lower flux and a lower 
SED peak energy than the Crab SSC emission, because the SED peak energy of the Vela 
pair SR is several decades lower. This difference is due to a combination of lower pair energy 
and lower Rlc- However, the Fermi spectrum can be produced by primary CR and SR. 

Simulating emission from two of the most energetic MSPs, B1821-24 and B1937+21, we 
find that their observed X-ray emission is well reproduced by SR from pairs. The SEDs of 
MSP SR components are predicted to peak at 1 - 10 MeV and could be detected by proposed 
Compton telescopes that would test this model. Since the soft X-ray emission is produced 
by the pairs at the low-energy turnover in the pair spectrum, well below the SED peak, the 
spectral index is that expected for SR from a single particle, —2/3 (4/3 for the SED). This 
index is much higher than the soft X-ray index of the pair SR spectrum of the Crab, where 
the soft X-ray emission is near the SED peak and produced by a range of pair energies. 
However, the predicted SSC fluxes of the MSPs are too low to be detectable with current 
instruments. 

Thus from our simulations, only Crab and Crab-like pulsars, such as B 1509-58 and 
B0540-69 in the LMC, are expected to have high enough levels of SSC emission at VHE 
energies to possibly be detectable by ground-based Air-Cherenkov telescopes. They share the 
characteristics that are important for production of strong SSC emission: high levels of non- 
thermal X rays, high pair multiplicity and high magnetic fields near the light cylinder. Non- 
Crab-like, middle-aged pulsars like Vela have much lower levels of non-thermal X rays relative 
to their GeV emission. So even if they are able to generate very high pair multiplicities ~ 10 5 , 
their pair energies will be too low and their SSC emission will not reach detectable levels 
at VHE energy. MSPs have higher levels of SSC than middle-aged pulsars, and their SSC 
spectra extend to higher energies than Crab-like pulsars, since their pair spectra reach TeV 
energies. However, their pair multiplicities may not be high enough to generate fluxes of 
SSC emission that are detectable by current telescopes. These findings are consistent with 
the recent stacking analysis of McCann (2015). 

Lyutikov (2013) presented a model for the Crab pulsar emission, assuming cyclotron 



and cyclotron-self Compton emission by pairs produced in the OG and gyrating at the 
cyclotron resonance with small pitch angles. It was suggested that the pairs acquire their 
pitch angles through coherent emission of radio waves, but their perpendicular momenta 
are assumed to remain non-relativistic so that the radiation occurs in the cyclotron rather 
than the synchrotron regime. The required pair multiplicity is M + = 10 6 — 10', well above 
what can be produced in either PC or OG pair cascades. In our model, the pairs gain pitch 
angles through resonant absorption and acquire relativistic perpendicular momenta, so that 
the radiation is in the synchrotron regime. Since the synchrotron radiation power is much 
higher than cyclotron power, we can produce the Crab optical to X-ray emission with a much 
lower pair multiplicity of 3 x 10 5 . 
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Fig. 1. — Pair spectra (in units of pairs/s/mc 2 ) from polar cap pair cascade simulations for 
the Crab, Vela, B1821-24 and B1937+21 pulsars, as labeled. The dashed line is a power law 
extension to the Crab pair spectrum (see text). 
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Fig. 2. — Evolution of the dynamics and radiation of pairs at initial energy y + as a func- 
tion of radius (in units of light cylinder radius) along their trajectory, for the Crab, Vela 
and B1937+21. Quantities plotted are perpendicular momentum, p± (in units of me), syn- 
chrotron loss rate, (d^//di) SR (cm -1 ), SSC loss rate, (d , y/d(!) ssc (cm -1 ), cyclotron absorption 
rate, (dp±/ dt) abs (cm -1 ), and critical synchrotron energy e gR (in units of me 2 ). 
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Fig. 3. — Evolution of the dynamics and radiation of primary electrons as a function of radius 
(in units of light cylinder radius) along their trajectory, for the Crab, Vela and B1937+21. 
Quantities plotted are Lorentz factor, 7 , synchrotron loss rate, (d'y/di) SR (cm -1 ), CR loss 
rate, (d^/d£) CR (cm -1 ), SSC loss rate, (d , y/d£) ssc (cm -1 ), and acceleration gain rate, i? acc 
(cm -1 ). 
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Fig. 4. — Model spectra of phase-averaged pulsed emission components from primary elec- 
trons and pairs (as labeled) from the Crab pulsar, for magnetic inclination angle a = 45° 
and observer angle ( = 60° and pair multiplicity M + = 2 x 10 4 . Data points are from Kuiper 
et al. (2001) [http://www.sron. nl/divisions/hea/kuiper/data.html|, Abdo et al. (2013) 
[http://fermi.gsfc.nasa.gov/ssc/data/access/lat/2nd_PSR_catalog/], Aleksic et al. (2012) 
and Alin et al. (2011). 
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Fig. 5. — Model spectra of phase-averaged pulsed emission components from primary elec- 
trons and pairs (as labeled) from the Crab pulsar, for magnetic inclination angle a = 45° 
and observer angle Q = 60° and pair multiplicity M + = 3 x 10 5 . The dashed lines are the 
SR and SSC spectra resulting from a power law extension to the cascade pair spectrum. 
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Fig. 6. — Model spectra of phase-averaged pulsed emission components from primary 
electrons and pairs (as labeled) from the Vela pulsar, for magnetic inclination angle 
a = 75° and observer angle £ = 60°. Solid lines are for pair multiplicity M + = 
6 x 10 3 and dashed lines are for M + = 10 5 . Data points are from Abdo et ah 
(2013) |http://fermi. gsfc.nasa.gov/ssc/data/access/lat/2nd_PSR_catalog/ and Harding et 
al. (2002). 
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Fig. 7. — Model spectra of phase-averaged pulsed emission components from primary elec- 
trons and pairs (as labeled) from PSR B1821-24, for magnetic inclination angle a = 45° and 
observer angle ( = 80°. Solid lines are for pair multiplicity M + = 10 3 and dashed lines are 
for M + = 10 5 . Data points are from Kuiper et al. (2005) and Johnson et al. (2014). The 
thick pink and yellow lines are the HESS and CTA sensitivity limits, respectively. 
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Fig. 8. — Model spectra of phase-averaged pulsed emission components from primary elec- 
trons and pairs (as labeled) from PSR B1937+21, for magnetic inclination angle a = 75° 
and observer angle £ = 70°. Solid lines are for pair multiplicity M + = 10 5 and dashed lines 
are for M + = 10 3 . Data points are from Ng et al. (2014) and Guillemot et al. (2012). The 
thick yellow line is the CTA sensitivity limit. 


